#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Merge huc12 and county to create data linkage

#load libraries
library(ggplot2)
library(raster)
library(sf)
library(dplyr)
library(readr)
library(mapview)
library(exactextractr)
library(data.table)
library(tidyverse)
library(Hmisc)
library(haven)
library(stargazer)
library(purrr)
library(tidyr)
library(collapse)
library(ggthemes)
library(mapview)
options(scipen=999)
library(stringr)

#open huc 12 data. only need MARB
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

MARB<-st_read("Data/WatershedBoundaries/Processed/MARB/MARB.shp")

#open county sf
#downloaded from https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
csf<-st_read("Data/CountyBoundaries/Raw/cb_2018_us_county_500k/cb_2018_us_county_500k.shp")


#make sure same crs
crs(MARB)
crs(csf)

#change county crs
csf<-st_transform(csf, crs=st_crs(MARB))
crs(csf)

table(csf$STATEFP)
table(csf$COUNTYFP)
csf$fips<-paste0(csf$STATEFP, csf$COUNTYFP, sep="")

csf<-subset(csf, select=c(fips, NAME, geometry))
names(csf)<-c("fips", "county", "geometry")

csf$fipsarea<-st_area(csf) %>% units::set_units("acres")

MARB$hucarea<-st_area(MARB) %>% units::set_units("acres")

#which hucs go to which fips?
link<-st_intersection(MARB, csf)

#find area
link$intarea<-st_area(link) %>% units::set_units("acres")

#create key
key<-subset(as.data.frame(link), select=c(huc12, fips, county, intarea, hucarea, fipsarea))
key<-key[order(key$huc12, -key$intarea),]
key<- key %>% group_by(huc12) %>% mutate(countn=n())
table(key$countn)

#9,412 huc12s are fully within one county. 
#less than half!

#proportion of huc 12 in county
key$prop<-as.numeric(key$intarea / key$fipsarea)
summary(key)

##will create weighted average with animal units data
#i will match animal units to the fips
#then, i will multiply animal units by proportion (huc 12 area/fips area)
#then, i will add the animal units by huc12. (collapse)

county_huc12_key<-key
save(county_huc12_key, file="Data/CountyBoundaries/Processed/Huc12County.RData")
